/****** generic program for the exogenous treatment case *******/
capture program drop test_exo
	program define test_exo, eclass 
	tempname diff 
	// testing a set of quantiles
	if testmean==0 {
	tempvar qtindicator
	matrix `diff'=J(1,numqt*(numl-1),.)
	local qind=0
	gen `qtindicator'=0
    foreach indqt of global qt {
	local qind=`qind'+1
	local indct=`indqt'*100
	qui centile y if t==1, centile(`indct')
	qui replace `qtindicator'=(y<=r(c_1)) if t==1
	qui centile y if t==0, centile(`indct')
	qui replace `qtindicator'=(y<=r(c_1)) if t==0
	local lind=0
		foreach l of global level {
		local lind=`lind'+1
		qui sum `qtindicator' if t==1 & x==float(`l')
		scalar temp=r(mean)
	    qui sum `qtindicator' if t==0 & x==float(`l')
		matrix `diff'[1,numqt*(`lind'-1)+`qind']=r(mean)-temp
	}
	}
	}
	// the mean test
	else if testmean==1{
	tempvar rank rank1
	matrix `diff'=J(1,(numl-1),.)
	cumul y if t==0, generate(`rank')
	cumul y if t==1, generate(`rank1')
	replace `rank'=`rank1' if t==1
		foreach l of global level {
		local lind=`lind'+1
		qui sum `rank' if t==1 & x==float(`l')
		scalar temp=r(mean)
	    qui sum `rank' if t==0 & x==float(`l')
		matrix `diff'[1,`lind']=r(mean)-temp
	}
	}
	ereturn post `diff'
	end 
	

/****** simulation set-up for the exogenous treatment case (Table A.1, DGP 1) *******/
/****** returns: bsstat (test statistic, where the variance covariance matrix is bootstrapped); */
/******          decchi (decison of the test, using the Chi2 critical value) *******/
/******          decf (decison of the test, using the F critical value) *******/
/****** decchi and decf will the same if the sample size is large *****/
capture program drop mcmain_exo 
program define mcmain_exo, rclass
	syntax [, obs(integer 1000) bsreps(real 500) testmean(integer 0) b(real 1)]
	set obs `obs'
	tempvar maxx 
	tempname dm wald hatv 
	gen x=2*ceil(runiform()*numl)/numl
	levelsof x, local(level)
	egen `maxx'=max(x)
	local maxx=`maxx'
	global level: list level-maxx
	gen t=round(runiform())
	gen v=rnormal()
	gen u0=rnormal()
	gen u1=rnormal()
	gen y1=x+1-`b'*v*x+v+u1 
	gen y0=x+v+u0	
	gen y=y0*(1-t)+y1*t
	scalar testmean=`testmean'
	// call the test function
	qui test_exo
	matrix `dm'=e(b)
	// call bootstrap using the test function to get the variance covariance matrix
	bootstrap _b, reps(`bsreps'): test_exo 
	matrix `hatv'=e(V)
	if testmean==0{
	scalar df=numqt*(numl-1)-diag0cnt(`hatv')
	}
	else if testmean==1{
	scalar df=(numl-1)-diag0cnt(`hatv')
	}
	return scalar df=df
	matrix `wald'=`dm'*invsym(`hatv')*`dm''
	return scalar bsstat=`wald'[1,1]
	return scalar decchi=(`wald'[1,1]>invchi2(df,0.95))
	return scalar decf=(`wald'[1,1]>invF(df,`obs',0.95)*df)
	drop _all
	end
	
/****** generic program for the endogenous treatment case *******/
	capture program drop test_end
	program define test_end, eclass 
	// testing a set of quantiles
	if testmean==0{
	tempname diff
	tempvar qtindicator
	matrix `diff'=J(1,numqt*(numl-1),.)
	tempvar y2
	qui gen `y2'=y
	qui replace `y2'=0 if t==1
	qui ivqte y (t=z), quantiles($qt) 
	matrix qte=e(b)
	qui ivqte `y2' (t=z), quantiles($qt)    // use y2 to trick ivqte estimate quantile under control
	matrix qt_c=-e(b) // notice the negative sign
	matrix qt_t=-e(b)+qte
	local qind=0
	gen `qtindicator'=0
    foreach indqt of global qt {
	local qind=`qind'+1
	qui replace `qtindicator'=(y<=(qt_c[1,`qind'])*(1-t)+qt_t[1,`qind']*t) 
	local lind=0
	foreach l of global level {
		local lind=`lind'+1
		qui sum `qtindicator' if z==1 & x==float(`l')
		scalar temp=r(mean)
	    qui sum `qtindicator' if z==0 & x==float(`l')
		matrix `diff'[1,numqt*(`lind'-1)+`qind']= r(mean)-temp
		}
	}
	}
	// the mean test
	else if testmean==1{
		tempname diff
		matrix `diff'=J(1,(numl-1),.)
		tempvar y2 rank
		qui gen `y2'=y
		qui replace `y2'=0 if t==1
		qui ivqte y (t=z), quantiles(0.01(0.01)0.99) 
		matrix qte=e(b)
		qui ivqte `y2' (t=z), quantiles(0.01(0.01)0.99)    // use y2 to trick ivqte estimate quantile under control
		matrix qt_c=-e(b) // notice the negative sign
		matrix qt_t=-e(b)+qte
		gen `rank'=0.5
		forv q = 1(1)99{
		qui	replace `rank'=`q'+0.5 if y>= (qt_c[1,`q']*(1-t)+qt_t[1,`q']*t)
		}
		foreach l of global level {
			local lind=`lind'+1
			qui sum `rank' if z==0 & x==float(`l')
			scalar temp=r(mean)
			qui sum `rank' if z==1 & x==float(`l')
			matrix `diff'[1,`lind']=r(mean)-temp
	}
	}
	ereturn post `diff'
	end 

/****** simulation set-up for the endogenous treatment case (Table A.1, DGP 2) *******/
/****** returns: bsstat (test statistic, where the variance covariance matrix is bootstrapped); */
/******          decchi (decison of the test, using the Chi2 critical value) *******/
/******          decf (decison of the test, using the F critical value) *******/
/****** decchi and decf will the same if the sample size is large *****/
capture program drop mcmain_end 
program define mcmain_end, rclass
	syntax [, obs(integer 1000) bsreps(real 500) testmean(integer 0) b(real 1)]
	set obs `obs'
	tempvar maxx 
	tempname dm wald hatv
	gen x=2*ceil(runiform()*numl)/numl
	levelsof x, local(level)
	egen `maxx'=max(x)
	local maxx=`maxx'
	global level: list level-maxx
	gen z=round(runiform())
	gen v=rnormal()
	gen u0=rnormal()
	gen u1=rnormal()
	gen y1=x+1-`b'*v*x+v+u1 
	gen y0=x+v+u0
	gen t=(0.15*(y1-y0)+z-0.5>0) 
	gen y=y0*(1-t)+y1*t	
	scalar testmean=`testmean'
	// call the test function
	qui test_end
	matrix `dm'=e(b)
	// call bootstrap using the test function to get the variance covariance matrix
	bootstrap _b, reps(`bsreps'): test_end
	matrix `hatv'=e(V)
	if testmean==0{
	scalar df=numqt*(numl-1)-diag0cnt(`hatv')
	}
	else if testmean==1{
	scalar df=(numl-1)-diag0cnt(`hatv')
	}
	return scalar df=df
	matrix `wald'=`dm'*invsym(`hatv')*`dm''
	return scalar bsstat=`wald'[1,1]
	return scalar decchi=(`wald'[1,1]>invchi2(df,0.95))
	return scalar decf=(`wald'[1,1]>invF(df,`obs',0.95)*df)
	drop _all
	end


